Light and stable triplet bipolarons on square and triangular lattices 
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We compute the properties of singlet and triplet bipolarons on two-dimensional lattices using 
the continuous time quantum Monte Carlo algorithm. Properties of the bipolaron including the 
total energy, inverse mass, bipolaron radius and number of phonons associated with the bipolaron 
demonstrate the qualitative difference between models of electron phonon interaction with long- 
range interaction (screened Frohlich) and those with purely local (Holstein) interaction. A major 
result of our survey of the parameter space is the existence of extra-light hybrid singlet bipolarons 
consisting of an on-site and an off-site component on both square and triangular lattices. We also 
compute triplet properties of the bipolarons and the pair dispersion. For pair momenta on the 
edge of the Brillouin zone of the triangular lattice, we find that triplet states are more stable than 
singlets. 



PACS numbers: 71.38.Mx, 71.38.Fp, 71.10.Fd 



I. INTRODUCTION 

Layered superconductors such as the cupratesi, iron 
pnictides^ and others such as hydrated cobaltates 
(Naa;Co02 • I.3H2O) (Ref. i) have the potential to 
completely change thinking about the origins of su- 
perconductivity. Given the importance of such a re- 
think, a number of competing theories and mechanisms 
have been suggested, which aim to understand supercon- 
ductivity in layered materials on at least a qualitative 
level^i^iii^i^iiSiiii. For example, it has been suggested 
that one way of explaining high temperature supercon- 
ductivity in cuprates is the condensation of pre-formed 
bipolarons^. A complete quantitative understanding of 
bipolaron formation, especially in the intermediate cou- 
pling regime is needed to judge the applicability of such a 
theory, especially when the electron-phonon interaction 
is not perfectly local (i.e. when the Holstein interaction 
is not appropriate) as is to be expected in real materials. 

There are a number of exact numerical methods 
for studying bipolaron formation in one-dimensional 
(ID) systems with electron-phonon interactions such as 
exact diagonalization (ED) ^^i^^ , advanced variational 
techniquesii, density matrix renormalization group 
(DMRG>1^ and various quantum Monte-Carlo (QMC) 
approachesiSiiiii^. In two-dimensional (2D) systems, 
many of those techniques are not applicable. Exact di- 
agonalization cannot cope with large numbers of lattice 
sites, DMRG schemes are not easy to develop in 2D and 
both ED and DMRG can suffer from a truncation of 
the phonon Hilbert space, rendering those methods in- 
efficient at strong electron-phonon coupling. Also, ad- 
vanced variational techniques are not easily generalized 
to 2D. On the other hand, QMC techniques (and espe- 
cially continuous time QMC) can cope with large lattice 
sizes and treat the phonon degrees of freedom exactly, 
making them well suited for computations of bipolarons 



on two- and potentially even three- dimensional lattices. 

There is a fascinating possibility relating to bipolarons 
in 2D. If a lattice is made up from triangular plaquettes, 
and if there is a strong Coulomb repulsion keeping bipo- 
larons from pairing on-site, superlight small bipolarons 
can form^^i ^i^^i^^ . Our aim here is to comprehensively 
compute the properties of bipolarons with long range in- 
teraction on both square and triangular lattices to de- 
termine if there are any other promising ways that bipo- 
larons can pair. The Hubbard-Holstein bipolaron on the 
square lattice was studied by Macridin et al. using the 
related diagrammatic quantum Monte Carlo technique 
(DMC)^'^. Our current study goes beyond this previous 
work by computing the properties of bipolarons formed 
from long range electron-phonon interactions. We also 
compute properties of bipolarons on triangular lattices, 
which have not been comprehensively studied even for 
the Holstein interaction. Our continuous-time quantum 
Monte Carlo (CTQMC) algorithm is particularly efficient 
for this task, since phonon degrees of freedom can be 
treated exactly to generate an effective retarded interac- 
tion between electrons. 

In this article, we study the screened Hubbard-Frohlich 
model, which has the Hamiltonian, 
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Here, c]^ creates an electron on site n, t is the inter- 
site hopping integral, U is the Hubbard repulsion, M is 
the ion mass, uj the ion oscillation frequency, C the ion 
displacement and P the ion momentum. The electron- 
ion force function is /, and has the form fm{n) = 

K [(m — n)^ + 1] ''^^ exp (— |m — n|/i?sc)- The screen- 
ing radius Rsc controls the length of the interaction, and 
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K is the strength of the interaction, m represents the 
position of the ion and n is the position of the electron. 
The electron-ion interaction leads to an effective retarded 
electron-electron interaction characterized by the func- 
tion, 

<i>A.(r(r),r(r')) = J2fm[rir)]frr.+Ar[r{r')] (2) 

m 

where Ar is the offset between end configurations of the 
paths. By changing Rsc it is possible to investigate the 
Hubbard-Holstein model {Rsc 0) and the Hubbard- 
Frohlich model {Rsc oo). We also consider a model 
similar to the nearest-neighbor model of Bonca and 
Trugman^i, where effective electron-electron interactions 
are truncated at near-neighbor sites ($(a) — $(0)/2 and 
all other $ = where a are vectors to the near-neighbor 
sites). This simplified model maps directly onto aU — V 
model in the large phonon frequency limit, and could be 
particularly useful for understanding limiting behavior—. 

This article is structured as follows. In section [III we 
briefly review the differences between the CTQMC algo- 
rithm in ID and 2D. In sections IlIII and [TVl we make a 
comprehensive survey of the parameter space of singlet 
bipolarons in 2D. Triplet properties are discussed in sec- 
tion|Vl One of the additional advantages of our algorithm 
is that bipolaron dispersions can be computed efhciently, 
especially in the large A regime, and these are shown in 
section IVTl Finally we summarize in section fVIII 

II. METHOD 

We use a continuous time quantum Monte Carlo al- 
gorithm based on path integrals. We have previously 
discussed our algorithm in detail with regard to compu- 
tations in ID^^, so we do not repeat those details here. 
There are a few small subtleties relating to using the al- 
gorithm in 2D, especially on a triangular lattice, which 
are discussed here. Our algorithm has been thoroughly 
checked against results from the U-V model. 

A. Triple kink insertions on triangular lattices 



kinks. Thus, to ensure ergodicity, an update with a 3 
kink insertion is proposed. 

Our scheme for 3 kink insertion is similar to the one 
that we use for binary updates. To avoid complications, 
we do not weight the positions of the kinks in imaginary 
time. Ternary kink insertions do not need to be especially 
efficient; if path configurations can be updated from even 
to odd numbers of kinks (and vice- versa) with reasonable 
regularity, binary insertions can be used to sample the 
remaining configurations efficiently. Our scheme is as 
follows, 

1. We select a kink type from the 6 possible kinks and 
assign the label h. 

2. We choose kinks at 120° and 240° rotations from 
li and assign them the labels I2 and ^3 

3. We choose insertion or removal of kinks with equal 
probability 1/2. 

4. If inserting, we choose imaginary times ri , T2 and 
for the new kinks with equal probability 1//3 from 
the interval [1, /?). 

5. If removal is selected and there is not at least one 
of each type of kink, then abort. Otherwise select 
a kink of type li for removal with equal probability 
1/Ni^, etc 

If configuration [D] has three more kinks than config- 
uration [C], then insertion takes place with probability 

P[C -^D]= min ( r^^^, ^e^[^l-^[^l, l) 

(3) 

and for removal, 

PiD ^ C] ^ mm(^'-[^]^'-g^'3[^]^.t.]-.[ci A 

[ [tpY J 

(4) 

Note that the number Ni-^[D\ represents the number of 
kinks of type ii in configuration D. Therefore when kinks 
are inserted, this is the number of kinks in the final con- 
figuration. When kinks are removed, Ni-^ [D\ is the num- 
ber of kinks in the initial configuration. 



There is a subtle difference to the algorithm on tri- 
angular lattices, or more generally on any lattice where 
an electron can return to its original position in 3 hops 
(such as lattices with nearest and next-nearest neighbor 
hopping). We illustrate this by considering an example 
configuration that is permitted on such lattices: Let one 
of the paths have no kinks and the other path have 3 
kinks - one in each of the nearest neighbor directions - so 
that the start and end of each path lies on the same lat- 
tice site. Such a configuration is clearly permitted, since 
the periodic boundary conditions in imaginary time are 
satisfied. However, if only binary kink insertions are in- 
cluded in the algorithm then it is not possible to update 
between configurations with odd and even numbers of 



B. Path exchange in 2D 

There is a small difference between exchange in ID and 
2D, especially on the triangular lattice. In ID, exchange 
can be carried out by inserting and removing kinks and 
antikinks. The form of the exchange update in ID is, 

P{C -^D) = min{l, Q[%Q^^^^i exp(A[C] - A[D])} (5) 
where, 

q(A) ^ ^.^.2nA-A niin(iVA-;,A) + l n^^.Pa-u^ 
mm{NAi + UA, A) + 1 

(6) 



3 



here A is the displacement between the t = (3 ends of 
path A and path B. ^Pk = n\ / (n — k)\ is the number of 
permutations. There are Np — min{NA-i, + l possible 
updates that can be made by inserting ua kinks into path 
A and removing niA antikinks, where ua + mA — A. ttia 
is chosen with equal weighting l/Np. Nai is the number 
of kinks of type I on path A. The expression for path 
B is similar, but the kink and anti-kink assignments are 
reversed to get displacement —A. 

In 2D, two sets of kinks and antikinks in different di- 
rections are required to exchange the ends of the paths, 
so the exchange update has the form, 

d 

P{C^D)=mm{^,exp{A[C]~A[D])l[Q[^l^^Qi%^} 

i=l 

(7) 

where there is a different A for each kink type to be 
inserted (and thus a different Np and m which should 
be chosen independently for each kink type and path) 
as determined from the number of kinks that would be 
required to hop from the end of one path to the end of 
the other, d is the number of lattice dimensions. 

For the triangular lattice, the choice of the two kink 
directions is not unique, since the kinks are not orthogo- 
nal. There are six nearest neighbor vectors representing 
kinks. To carry out the update, kinks representing two 
different directions are chosen with equal weight from all 
available possibilities (that is the two chosen kinks may 
not be antikinks of each other). Since an equal weight 
scheme is used, the factors relating to the probability of 
this choice cancel on both sides of the balance equations. 
It is important to choose kinks from all possible direc- 
tions to improve the efhciency of the algorithm. 

C. Global path shifts 

As in the ID case, a useful update displaces a single 
path when the configurations are not exchanged. In both 
cases, one of the paths is chosen with equal probability, a 
kink direction is chosen with equal probability and then 
the path is shifted through r lattice sites, where r is an 
integer chosen randomly between and i?max ~ 1- Typ- 
ically i?inax = 50 for lattice size 100. Failure to include 
global path shifts leads to an inefficient algorithm with 
large correlation between measurement, which is espe- 
cially bad when the bipolaron is only just bound. Global 
shifts are essential when the bipolaron radius is to be 
measured. 



III. SINGLETS ON THE SQUARE LATTICE 

The survey of the bipolaron parameter space be- 
gins by considering singlet properties of bipolarons on 
the square lattice. We will examine the effects of the 
electron-phonon coupling, Hubbard repulsion and inter- 
action range on the bipolarons. In the following, we 



take uj = Lojt = 1, which corresponds to ujjW = 1/4 
which is well inside the adiabatic regime. W is the half 
band width or energy of a single non-interacting electron. 
W — At for the square lattice and W = 6t for the trian- 
gular lattice. We restrict paths to lie within 100 lattice 
spacings of each other, with the large lattice size used to 
ensure that the effects of the lattice restriction are min- 
imal (this should be true where the bipolaron radius is 
much smaller than the lattice size). Most properties are 
not strongly affected by finite size effects. We note that 
the inverse radius is the most sensitive measurement, and 
the total energy the least sensitive. The most severe fi- 
nite size effect is that premature binding can occur if the 
paths are too closely restricted. 

Figure [1] shows the total energy of the two polarons 
/ bipolaron pair for a range of electron-phonon coupling 
strengths, interaction screening radius and Hubbard U. 
As is the case in IDi^, a qualitative difference can be seen 
between bipolarons formed with Holstein and screened 
Frohlich interactions. At large U, the energy curves for 
the Hubbard-Holstcin bipolaron arc flat, indicating that 
the bipolaron has unbound into two polarons. This is 
true for all the A shown, indicating that there is no bound 
singlet bipolaron in the [/ — > oo Hubbard-Holstein model. 
The significance for triplet states will be discussed later 
in the article. 

In contrast to the flat energy curves characterizing 
the large U, large A Holstein model, the energy of the 
Hubbard- Frohlich bipolaron changes on varying U , indi- 
cating that there are bound singlet states even for very 
large U (we will revisit this point later in the article). 
For small A and Rsc ^ 0, the energy curves are flat at 
large U. The model with only nearest-neighbor interac- 
tion {^{a) — 0.5) has similar properties to the Rsc — 2 
and Rsc — 3 models. $(a) = 0.438 to 3 significant fig- 
ures for the model with Rsc = 2 and <f>(a) = 0.509 to 3 
significant figures for the model with Rsc — 3, indicating 
that bipolaron properties depend most strongly on the 
nearest neighbor part of the Frohlich interaction, rather 
than the tails. 

Is is useful to quickly describe the limiting behaviors 
of the bipolarons. First, it is appropriate to introduce 
some notation. We use SO to define a strongly bound 
on-site singlet pair where the two electrons are separated 
by no lattice spacings and use SI to denote a strongly 
bound inter-site pair where the electrons are separated 
by one lattice spaced. As discussed for the ID modcl^-, 
if bipolarons are strongly bound the renormalizcd inter- 
site hopping is small, and an atomic Hamiltonian can be 
determined using the Lang-Firsov transformation: 

Tin' \ 1 / m ^ ' 

(8) 

Since the electron number operator, n is unchanged un- 
der the Lang-Firsov transformation, the effects of the 
Hubbard U can quickly be reintroduced. Thus, for a 
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FIG. 1: Total energy of bipolarons on the square lattice for various U, X, Rsc, u) = 1 and /3 = l3/t = 14. There is a 
clear qualitative difference between bipolarons with long range interactions and those formed from the purely local Holstein 
interaction. Unless otherwise indicated, error bars represent three standard deviations. Where error bars are not visible, they 
are smaller than the points. 



0.4 
0.3 
0.2 
0.1 






Square, Frohlich 












□ H 


g □ □ □ □ □ 


□ □ O □ □ E 


!! 4- • 








0.5 
0.4 
0.3 

0.2 
0.1 


0.5 
0.4 
0.3 
0.2 
0.1 



Square, Frohlich, Rsc=3 










„-^3--B- H.--H □ □ □ 


D □ E E 







0.5 1 1.5 2 2.5 3 3.5 

UIW 



jm [uniiinDaQE a-B □HEDnHmciiH- 



Square, Holstein 



a-e-e-9© e © e .? 




Square, Frohlich, R^q=2 


-X— -- --X— --X 








X' 








-.-"^ _ H s Q ° ^ ° 


■ □ ■ □ □ -E 








FIG. 2: Inverse mass of bipolarons on the square lattice. Parameters are the same as in Fig. [T] Again, qualitative differences 
between the bipolarons formed via the purely site local Holstein interaction and the long-range interactions can be seen. 



strongly bound onsite (SO) bipolaron, 

E^u~ m\ (9) 

W is the non-interacting kinetic energy of a single par- 
ticle (half band- width) . The limiting behavior of the SO 
bipolaron is shown on the plots as dashed lines. As in 
the ID case, the SO Hubbard-Holstein bipolaron rapidly 
becomes strongly bound on decreasing \J . For the longer 
range interactions, it takes much more attraction to 
strongly bind the SO bipolaron. We can also compute the 



energy of the strongly bound SI bipolaron which forms 
at large A and U » W\, 

i; = -2VFA[l + $(a)/$o]- (10) 

The energy of the strongly bound SI state is plotted in 
Fig. [1] as the arrows at U /W — 3.5 (note that only 5 
arrows corresponding to the largest A are shown). The 
energy of the strongly bound SI bipolaron is higher than 
the simulated energy, but the values converge as A in- 
creases. From the difference in energies between the SI 
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FIG. 3: Comparison of the inverse mass for bipolarons on 
the square lattice. A = 1.6 and various interaction ranges are 
considered. Other parameters are as in Fig. [1] At interme- 
diate U, the energies of the on-site (SO) and intersite (SI) 
bipolarons are degenerate, so a light bipolaron can be formed 
without the need for triangular plaquettes. 



limit and the simulated energy, it appears that the bipo- 
larons formed from heavily screened interactions have a 
stronger SI character. 

The inverse mass of the bipolarons is shown in Fig. [21 
Again, significant differences between Hubbard-Holstein 
and Hubbard-Frohlich bipolarons can immediately be 
seen. In the Hubbard-Holstein model, the mass changes 
very rapidly with U when the bipolaron binds, with a 
mass change of several orders of magnitude over a very 
small range of U values. This rapid change is not seen for 
the Frohlich bipolaron, which maintains a similar mass 
over a wide range of U. For the larger A values and large 
U, the system of two particles appears to be heavier in 
the Frohlich case. This is because the Holstein bipolaron 
is not bound for those values, whereas the Frolich inter- 
action leads to an intersite (SI) bipolaron. 

In comparison to the U dependence of the mass of the 
ID Holstein bipolaron, which can be found in Ref. [TgI . 
the change in mass with U is extremely rapid on the 
square lattice. In contrast, the mass relating to long 
range interactions is qualitatively unchanged, with the 
formation of a relatively light bipolaron (even in com- 
parison with a free electron) over a very wide range of 
the parameter space. The qualitative difference between 
the properties of Holstein and Frohlich bipolarons em- 
phasizes that while the Holstein model may be a reason- 
able approximation for considering polarons in organic 
molecular compounds, the simple Holstein approxima- 
tion is inappropriate for many crystalline systems where 
complete screening of the electron-phonon interaction is 
not possible. Even if the screening length is as small 
as a single lattice spacing, bipolarons have significantly 
modified properties compared with those formed in the 
Holstein model. 

In ID, a crab like bipolaron can form in models with 
long range interaction when the energies of on-site and 
off-site configurations are degenerate (SI like configura- 



tions can be changed to SO with a single hop, without 
any energy barrier to make the process second order in 
t) . This leads to a noticeable decrease in mass at interme- 
diate U. It is of interest to determine if these lighter bipo- 
larons can also form on the square lattice. In Fig. [3l the 
masses of bipolarons formed when A — 1.6 are examined 
as Rsc is changed. A peak in the inverse mass (decrease in 
the mass) is clearly visible at intermediate U and is espe- 
cially pronounced when Rsc = 2. Comparing equations [9] 
andfTUl it can be seen that for large phonon frequencies, 
the peak would be expected at U/W = 2A[l-$(a)/$(0)], 
which corresponds to U/W = 1.798 for Rsc = 2 (to 4 s.f) 
in good agreement. The presence of this peak shows that 
light bipolarons can exist on square lattices without the 
need for triangular plaquettes. Naturally, such a state de- 
pends on a very subtle balance of parameters and could 
not be realized in all crystalline systems. 

As in the case of the ID bipolaron, analytical deter- 
mination of the effective mass, even for limiting values 
of the parameters, is difficult (and may not be possible) 
when there are long range tails in the interaction. While 
we have not computed limiting behaviors of the mass for 
the model with only near-neighbor interaction, we ex- 
pect that a similar approach to that applied to the ID 
bipolaron by Bonca and Trugman could yield results^^. 

We also show the number of phonons associated with 
the bipolaron cloud (Fig. U]). As before, the properties 
of the Hubbard-Holstein bipolaron are qualitatively dif- 
ferent to the bipolarons formed from long-range electron 
phonon interaction. The number of phonons associated 
with the Holstein bipolaron shows an abrupt increase as 
the bipolaron binds on decreasing U, whereas there is 
a smoother crossover from SI to SO behavior associated 
with the screened Frohlich interaction. 

The total number of phonons associated with the 
strongly bound onsite bipolaron (following the argument 
in Ref. M) is. 
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(11) 



This number is independent of U, so the number of 
phonons reaches a limiting value on decreasing U as seen 
in Fig. m (the strongly bound SO limit is represented by 
arrows on the left of the graph). The Holstein bipolaron 
with large A quickly reaches the SO limit at low U. Bipo- 
larons formed with long range interactions approach the 
SO limit less rapidly since the longer range tails lead to a 
shallower effective potential. 

The number of phonons associated with the SI bipo- 
laron can also be found, 
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This is represented as arrows on the large U side of the 
plots. The agreement improves on increasing A. Again, it 
is clear that the most heavily screened interactions have 
the most clearly defined SI character. 
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FIG. 4: Number of phonons associated with bipolarons on the square lattice. Parameters are identical to Fig. [T] Again, the 
properties of the Holstein bipolaron are qualitatively different to those of the Frohlich bipolaron. 
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FIG. 5: Inverse bipolaron size associated with bipolarons on the square lattice. Parameters as in Fig. [T] 



We complete this section by examining the inverse 
bipolaron size associated with the square lattice. Again, 
qualitative differences exist betvi^een the Holstein bipo- 
laron and the Frohlich bipolarons, with SI bipolarons 
forming in the screened Frohlich model at large U and A 
indicated by a radius of around a lattice spacing. Com- 
parison with figure [2] shows that bipolarons on the square 
lattice with long range interaction may be both small 
and light at intermediate and large A and intermediate 
U, where polarons are boimd into the hybrid SO-Sl bipo- 
laron. 



IV. SINGLETS ON THE TRIANGULAR 
LATTICE 



In this section, attention is turned to the properties 
of bipolarons on the triangular lattice, which have al- 
ready been shown to have unusual propertiea^i. To en- 
sure that results for bipolarons on square and triangu- 
lar lattices can be meaningfully compared, we have kept 
energy scales fixed in terms of W (the non-interacting 
kinetic energy of a single particle) rather than t. For the 
square lattice, W = At, while for the triangular lattice, 
W = 6t. Similar considerations should be made when 
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FIG. 7; Number of phonons associated with bipolarons on the triangular lattice. Parameters are as in Fig. (6] 



comparing plots of total energy between lattice types, 
which is why we have quoted all energy scales in these sec- 
tions in terms of W. Also, to keep the ratio uj/W = 1/4, 
all properties in this section are computed for uj/t = 1.5. 
Finally, for the same reason, the temperature is also 3/2 
higher than that used for computations of the bipolaron 
on the square lattice, so ^ = 28/3. 

Figure [6] shows the total energy of bipolarons on the 
triangular lattice as U, X and R^c are changed. A quick 
inspection indicated that there is no clear qualitative dif- 
ference between the total energy of bipolarons on square 



and triangular lattices. The same is true of the num- 
ber of phonons associated with the bipolaron, which is 
shown in Fig. [T] It is quite surprising that the prop- 
erties are so similar, since there are very different hop- 
ping mechanisms on the triangular lattice compared to 
the square lattice (such as the possibility for SI bipo- 
larons to move with a single hop on the triangular lat- 
tice, while intersite bipolarons on the square lattice have 
to tunnel through a high energy intermediate state to 
move). However, careful comparison shows that in the 
strong U limit (where inter-site SI bipolarons are formed) 
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the bipolaron formed from near-neighbor interactions has 
slightly more phonons on the square lattice than on the 
triangular lattice. For the very long range Frohiich in- 
teraction, there are only tiny differences between results 
for the square and triangular lattices (there are around 
2.5% more phonons associated with Frohiich bipolarons 
on the triangular lattice) . Clearly, the crossover between 
SO and SI bipolarons occurs at approximately the same 
value of U /W. The strong similarities between the prop- 
erties of bipolarons on square and triangular lattices on 
the adiabatic side of the parameter space is quite surpris- 
ing considering that the properties are very different in 
the antiadiabatic limit. However, the differences between 
the bipolarons becomes clearer when the inverse mass is 
examined. 

It has previously been reported that bipolarons on 



the triangular lattice are small and light^i. Figure [5] 
shows the inverse mass of bipolarons on the triangular 
lattice. Bipolarons formed via the Frohiich interaction 
are relatively light, as was observed when the square 
lattice was considered. The Hoistein bipolaron remains 
qualitatively different to the bipolarons formed with the 
screened Frohiich interaction, with an enormous change 
in the mass of several orders of magnitude over a very 
small range of U/W as the bipolaron binds. The most 
significant (and surprising) differences between the bipo- 
larons on square and triangular lattices can be seen when 
nearest neighbor interactions are considered and when 
A = 0.6 and U is large. For these parameters, the bipo- 
laron on the square lattice is lighter than its counterpart 
on the triangular lattice (although this does not mean 
that it is necessarily small). We also note that as A ap- 
proaches 1.6, the relative magnitudes of the masses on 
square and triangular lattices reverses and the bipolaron 
on the triangular lattice becomes lighter than the bipo- 
laron on the square lattice (the other parameters are held 
fixed). This is the part of the parameter space where su- 
perlight small bipolarons are expected. We will revisit 
this point when the size of the bipolaron is computed. 

Again, a hybrid SO-Sl bipolaron is likely when the SO 
and SI configurations become degenerate. As before, this 
state is visible as a decrease in the effective mass (increase 
in inverse effective mass) at intermediate U. A compar- 
ison of the inverse masses of bipolarons as interaction 
range is changed is shown in Fig. [51 There is only a shal- 
low maximum in the inverse mass curve for Rgc — oo, 
but there are more pronounced humps in the curves for 
smaller values of the screening, corresponding to a light 
bipolaron. The peaks are not as pronounced for the bipo- 
laron on the triangular lattice as they were when motion 
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FIG. 10: Inverse size associated with bipolarons on the triangular lattice. Parameters are as in Fig. [6l 



was on the square lattice. This is probably because the 
SI bipolaron on the triangular lattice at large A and large 
U is already light and small with a mass that is 1st order 
in t^, so the presence of the hybrid state does not make 
it significantly easier for the bipolaron to move about on 
the triangular lattice (whereas on the square lattice, the 
change from SI to the hybrid state changes bipolaron 
hopping processes from 2nd to 1st order in t leading to a 
reduction in mass). 

To finish the examination of the differences between 
square and triangular lattices, we examine the inverse 
bipolaron size in figure [TOl For large A, the size of the 
pair is small and SI bipolarons are clearly formed in the 
Frohlich model at large U (i.e. the radius is approxi- 
mately one lattice spacing). There are differences be- 
tween the sizes of bipolarons on triangular and square 
lattices. Comparison between the radius of the bipolaron 
on the two lattices shows that the bipolaron on the trian- 
gular lattice is significantly smaller than its counterpart 
on the square lattice at large U. Also, the bipolaron on 
the triangular lattice is slightly bigger than the pair on 
the square lattice at very small U. The bipolarons formed 
on square and triangular lattices have similar size at in- 
termediate U where the pair has hybrid SO-Sl properties. 
Thus we can see that the slightly larger (but still light) 
mass of the SI bipolaron on the triangular lattice is most 
likely a result of the stronger binding into the SI state 
(indeed for large U, the SI state is very well bound on the 
triangular lattice at quite small A - shown by R^^ « 1). 
It is this property of small and light bipolarons that could 
lead to a bipolaron condensate of SI bipolarons on a tri- 
angular lattice at reasonably high temperatures. How- 
ever, we have now identified an additional bipolaron con- 
figuration which is light and small on both triangular and 
square lattices, and which is formed at moderate A and 



U: the extra-light hybrid SO-Sl bipolaron. 



V. TRIPLETS 

We now direct our attention to triplet properties of 
the bipolaron. The possibility of triplet superconductiv- 
ity has received a lot of interest since the discovery of 
spin triplet superconductivity in Sr2Ru04^^. Triplet su- 
perconductivity has also been identified in heavy Fermion 
materials that have a triangular lattice^!. From the the- 
oretical point of view, triplet pairs could be of interest 
for two reasons. First, there is a wide literature on BCS 
to BEC crossove r^^i^^'^" which is of interest when the pa- 
rameters governing pairing are intermediate. The pres- 
ence of stable real-space triplet pairs would add addi- 
tional limits to this problem, including a crossover or 
transition between singlet and triplet pairing. Second, 
the energy difference between singlet and triplet states 
could be considered as the energy cost of a spin flip, 
with a direct interpretation as a spin gap'^^. In the next 
two sections of this article, we discuss the possibility and 
properties of real space triplet pairs on square and trian- 
gular lattices. 



A. Hubbard-Holstein model 

We start our examination by considering the on-site 
Hubbard-Holstein interaction. To examine if triplet bipo- 
larons can exist in the 2D Hubbard-Holstein model, we 
examine the model with very large U/W = 10 and 
compare the energy of two electrons subjected to the 
Hubbard-Holstein interaction with that of two unbound 
polarons (Fig. \TT^ . There are several good reasons for 
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FIG. 11: Absence of triplet bipolaron in the Holstein model 
on square and triangular lattices, (a) Square lattice, /3 = 7, 
id = 1, U/t = 40. (b) Triangular lattice, /3 = 14/3, uj = 
1.5, U/t — 60. Properties are computed at various A. Note 
that with the energy expressed in terms of the band-width, 
lattice type has very little effect on Holstein polarons. The 
singlet and triplet curves lie under the polaron curve, and are 
identical to within the statistical error. 

doing this. First, as we have seen in the previous sections, 
the energy of the singlet bipolaron increases monotoni- 
cally with U, so the binding energy of the singlet bipo- 
laron is minimized, therefore the singlet triplet splitting is 
also minimized (since triplet states can not have a higher 
energy than two free polarons without dissociating) . Sec- 
ond, since triplet states must have equal or higher energy 
than singlets, an unbound singlet also implies that the 
triplet is also not bound. The triplet properties can not 
depend on U, since a node in the triplet pair wavefunc- 
tion means that the Hubbard U has no possible effect on 
triplet properties, so it is sufficient to analyze only this 
large U limit. We have carefully computed the energies 
of singlets, triplets and two polarons on triangular and 
square lattices, and have found no bound triplet states 
to within the statistical error. 



B. Hubbard-Frohlich model 

While there are no triplet bipolarons formed in the 
Hubbard-Holstein model, the longer range interactions 
associated with the Hubbard-Frohlich model can lead 
directly to inter-site pairing, which is likely to lead to 



triplet states. To determine the existence of triplet pairs 
in the model, we compute the singlet and triplet energies, 
and compare them with the energy of two polarons, which 
can be seen in Fig. [T^ The parameters associated with 
the computations are U — AOt, /3 — 7 and cZ; = 2 on the 
square lattice and U/t — 60, Cj — 3 and (3 = 14/3 on the 
triangular lattice. We use the same parameters through- 
out this section. At weak A, there are neither singlet nor 
triplet bound states. On increasing A, the singlet and 
triplet states bind at a similar value of A (the binding 
can be seen as bifurcations of the curves). At very large 
A, the energies of the singlet and triplet become degen- 
erate. This is a result of the interaction on the polaron 
hopping, which becomes exponentially small. The sepa- 
ration of singlets and triplets in the U — V model on the 
square lattice scales as t^, whereas on the triangular lat- 
tice the large A singlet-triplet splitting scales like t. This 
can be seen in the adiabatic limit as a much larger singlet 
triplet splitting in the case of the triangular lattice. We 
note that there are strictly two degenerate triplet states 
measured here, representing two p-pairings. 

The inverse mass of the triplet Frohlich bipolaron is 
shown in Fig. [T31 At this value of U /W and large 
electron-phonon coupling, the triplet bipolaron is heav- 
ier than the singlet bipolaron, and at weak coupling the 
triplet is lighter. As U is changed, the triplet properties 
remain constant, whereas the singlet ones change, no- 
tably with an increase in the singlet mass on degrease in 
U. At very small U, we would expect that the singlet is 
always heavier. While the strongly coupled singlet bipo- 
laron is clearly lighter on the triangular lattice, which is 
due to the qualitative effects of lattice type, we do not 
see any qualitative difference between the mass of the 
triplets on square and triangular lattices which are very 
similar. 

Figure [14] shows the number of phonons associated 
with the triplet Frohlich bipolaron. At strong coupling 
on the triangular lattice, there is a clear difference be- 
tween the number of phonons associated with singlet and 
triplet. In contrast, on the square lattice the number of 
phonons associated with the two types of bipolaron con- 
verges on a single value at large A. We believe that this 
is related to the second order hopping on square lattices 
at strong coupling and first order hopping on triangu- 
lar lattices. We also note that increasing the electron- 
phonon interaction range reduces this difference, presum- 
ably since the bipolaron has more configurational free- 
dom when Rsc is increased, thus allowing similar hopping 
processes. 

Finally, we compute the radius of the triplet Frohlich 
bipolaron, which is shown in Fig. 1151 There are only 
small quantitative changes between the radii of triplet 
bipolarons on different lattices. We note a small finite 
size effect at weak coupling, which corresponds to a fi- 
nite inverse radius where the bipolaron is not bound. On 
increasing A, the singlet bipolaron is the first to bind 
strongly, with the triplet radius decreasing slightly, but 
with the triplet becoming small only for larger A. In- 
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on both lattices. However, the singlet-triplet splitting on the triangular lattice is much larger than on the square lattice for all 
interaction ranges. 



creased interaction range decreases the electron-phonon VI. DISPERSION 

coupling required to bind both types of bipolaron. The 

strongly bound triplet is slightly larger than the singlet. a. Hubbard-Holstein model 



To complete our survey of bipolarons in 2D, we con- 
sider dispersions of the Hubbard-Holstein and Hubbard- 
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At low A, the triplet bipolaron is lighter than the singlet bipolaron. Note that since properties of the triplet bipolaron do not 
depend on U, then the singlet will be always be heavier at low U. 



Frohlich bipolarons. The electron-phonon interaction has 
the potential to lead to unusual effects. For example po- 
laron dispersions in the adiabatic regime of the Holstein 
model on square and triangular lattices are flattened at 
the edge of the Brillouin zone^^i^. It is not known if 
the dispersion remains flattened once bipolarons form, 
and also little is known about the dispersions of triplet 



bipolarons. Thus, it is of interest to determine how the 
flat dispersion evolves as a bipolaron is bound from two 
polarons. 

In figure [16] we show singlet dispersions of the 
Hubbard-Holstein bipolaron on the square lattice when 
A — 1.45 and w = 1. Computation of the dispersion is 
carried out using the same method as in Ref. At large 
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FIG. 14: Number of phonons associated with the triplet Frohhch bipolaron. Parameters are as Fig. 1121 



J7, the bipolaron is unbound and the dispersion repre- 
sents two polarons. The dispersions relating to U jt — 12, 
\J jt = 11.6 and I] jt — 11.2 are indistinguishable to the 
eye (top panel). As the repulsion is decreased, there is a 
decrease in bandwidth of several orders of magnitude at 
the V corresponding to the binding of the bipolaron. In 
the lower panel, we also show dispersions normalized by 
Em — Co (^M is the energy of the bipolaron at the M point) 
so that the forms of the dispersion function can be com- 



pared to the non-interacting tight-binding dispersion. At 
large [/, the form of the dispersion of the two polarons 
(and the weakly bound bipolaron) is quite distorted away 
from the form of the tight-binding dispersion. As C/ de- 
creases further, the electrons become tightly bound into 
an on-site bipolaron. Then the form of the dispersion 
tends towards that of the tight-binding spectrum. We 
believe this to be the effect of the strongly bound bipo- 
laron acting as a single particle, with a new hopping in- 
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FIG. 15: Radius of the Frohlich triplet bipolaron. Parameters are as Fig. 1121 



tegral that relates to the overlap of the wavefunction of 
the pair between sites. 

In figure [TTI we show the variation of the singlet disper- 
sions of the Hubbard-Holstein bipolaron on the triangular 
lattice when A — 1.45, u) — 1.5 and U is varied. Again, 
the functional form of the spectrum tends towards the 
tight-binding dispersion as pairs become strongly bound, 
although it becomes increasingly difficult to collect data 
since the fractional variance on the dispersion increases as 



A increases. Again, we suggest that the similarity in func- 
tional forms is due to the tightly bound bipolaron acting 
as a single particle with a single (albeit small) inter-site 
hopping parameter. For the Holstein bipolaron, the pair 
wavefunction is small, and the lattice distortion is highly 
localized. 
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FIG. 16: Singlet dispersions of the Hubbard-Holstein bipo- 
laron on the square lattice when A = 1.45, /3 = 14 and oj = 1. 
U is varied. Note that error bars on the plots of dispersions 
represent one standard deviation. The lower panel shows the 
dispersions normalized by em — eo compared to the normal- 
ized tight binding dispersion. As U decreases, the bipolaron 
becomes strongly bound on-site and the dispersion develops 
a form closer to the tight-binding approximation. 



FIG. 17: Singlet dispersions of the Hubbard-Holstein bipo- 
laron on the triangular lattice when A — 1.45, (3 — 28/3 and 
= 1.5. U is varied. Error bars represent one standard er- 
ror. There is a very rapid change in the bandwidth between 
U/t = 15.6 and U/t = 16.8. This coincides with a change 
in the shape of the dispersion from the flattened to the tight 
binding form. 



B. Hubbard-Frohlich model 

We also examine the dispersions of the Hubbard- 
Frohlich model. In Fig. [THl we show the dispersion 
for the square lattice when the Frohlich interaction with 
Rsc = 2 is considered. The panels show the effects of 
changing U. At large U, the bipolaron is bound into 
an intersite (SI) configuration. There are two main fea- 
tures that should be emphasized. The first is remarkable. 
As the Coulomb repulsion decreases, the bandwidth in- 
creases significantly, i.e. although the bipolaron becomes 
more strongly bound, it becomes lighter. This is the same 
effect as seen in section 1111) where the mass was shown 
to decrease at intermediate U as a, new type of hybrid 
SO- SI crab bipolaron forms. At small U, the degeneracy 
of the SO and SI bipolarons is broken and the bandwidth 
decreases with U. 

The second feature in Fig. [TH]is that the shape of the 
bipolaron dispersion changes dramatically as it it bound 
onsite. For large U, the dispersion is flattened, whereas 
for the strongly bound onsite pair {U/t — 5) the dis- 
persion has the characteristic tight-binding shape (albeit 



strongly renormalized by interactions). We consider this 
effect to relate to the change between hopping mecha- 
nisms. For a weakly bound bipolaron, the wavefunction 
is large and single polaron-like hops dominate the motion 
of the particle, whereas for the strongly bound bipolaron, 
the pair wavefunction is smaller than a single lattice site 
and a new hopping integral becomes relevant, which is 
relating to the tunneling of the whole bipolaron between 
sites. The dispersions of triplet states on the square lat- 
tice are also shown in the panel relating to U/t — 40. 
The triplet state is higher in energy and has a narrower 
bandwidth. Singlet and triplet dispersions do not cross. 

We also examine the dispersions of the Coulomb- 
Frohlich bipolaron on a triangular lattice, which are 
shown in Fig. 1191 The first significant difference between 
bipolarons on square and triangular lattices is the exis- 
tence of super light bipolarons at large U ^ oo. These 
contribute to a factor of ~ 4 difference between the band- 
widths of the bipolaron on square and triangular lat- 
tices. Again, as U decreases there is a remarkable in- 
crease in the bandwidth. This effect is not as dramatic 
as in the case of the square lattice, presumably because 
the bipolaron already moves with a crab like motion at 



16 



0.25 I 

(a) 

0.2 



S 0.15 
d 0.1 
0.05 - 




singlet 
triplet 

Square, U/t=40 tight-binding 





k[a'] 



k[a' 



0.007 
0.006 
0.005 
0.004 
0.003 
0.002 
0.001 




(c) . 




singlet — ■ — ' \ 




tigtit-binding ----- ^ \^ 


\ 




\ 


/ Square, £//(=5 


\ 




V 




FIG. 18: Singlet dispersions of the Hubbard- Frohlich model on the square lattice. The triplet dispersion is shown in panel (a). 
For comparison, the tight-binding dispersion normalized to the bandwidth is also shown. A = 1.45, w = 1, Rsa = 2 and /3 = 14. 
Error bars show one standard deviation. As U decreases and the bipolaron becomes strongly bound onsite, the shape of the 
dispersion becomes consistent with that of a tightly-bound particle. An increase in the band-width can be seen initially on 
decreasing U corresponding to the formation of the hybrid SO-Sl bipolaron. For smaller (7, the SO bipolaron is tightly bound 
and the dispersion is narrow. 



large U . However the degeneracy between S'O and S\ 
bipolarons accounts for an additional ~20% increase in 
the bandwidth, making the bipolarons extra-hght. This 
shows how bipolarons can be light over an extremely large 
regime of the parameter space on triangular lattices. 

To highlight the effect of the formation of the hybrid 
SO-Sl bipolaron on the dispersion. Fig. shows a com- 
parison of the variation of the bandwidth as V is changed. 
The bandwidths of bipolarons on both square and trian- 
gular lattices are shown, although it should be noted that 
the scales are different. There is an increase in bandwidth 
at U jW ~ 3. The widening of the bandwidth is more sig- 
nificant for the square lattice, but is also visible for the 
triangular lattice. The broadening of the bandwidth on 
the square lattice at intermediate to large U is important, 
since it shows that light bipolarons are not confined to 
lattices constructed from triangular plaquettes, but can 
also exist on simple square, linear and presumably cubic 
lattices. 

As noted in the introduction, we can also consider the 
properties of bipolarons formed from a simplified nearest- 
neighbor interaction as a means to compare with analyt- 
ical results. Fig. [^T] shows the variation of the singlet 
and triplet dispersions in the near-neighbor model as lo is 



changed . Bipolarons are considered on both square and 
triangular lattices with a large electron-phonon coupling 
of A = 8. To make comparisons with the t/V^- model, we 
choose a large phonon frequency lo = VlQt on the trian- 
gular lattice (and w — 80i on the square lattice). A large 
coulomb repulsion oi U = 12Qt [U — %Qt) stops on-site 
pairing, which leads to a, U — V model with U = 2At 
{U = 16t) and V = 48t {V = 32t). We also compute the 
dispersions of bipolarons on the triangular lattice when 
uj — 30 and uj — 6 {Co = 20 and uj = A on the square 
lattice). 

The dispersion of the near-neighbor bipolaron on the 
triangular lattice at large phonon frequency yields an un- 
usual result. The singlet and triplet bands cross. We note 
again that the QMC algorithm picks out only the low- 
est energies associated with singlet and triplet states. In 
the plot, the s band is visible, along with the lowest p 
band. A crossing is very unexpected, since in all situa- 
tions that we have previously studied, the triplet band 
has always sat above the singlet one. This indicates that 
there will be a transition from singlet to triplet states if 
sufficient momentum can be imparted to the bipolarons. 
The crossing remains down to lower phonon frequencies 
huj 5W, but is not visible in the adiabatic regime. 
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FIG. 19: Singlet dispersions of the bipolaron formed in the screened Hubbard-Frohlich model on the triangular lattice. Rsc ~ 2, 
A = 1.45, ijj = 1.5, /3 = 28/3 and U is varied. Error bars represent one standard deviation. As U is decreased, the bandwidth 
initially increases as the SO-Sl hybrid bipolaron is approached from the SI bipolaron (until around U/t = 15). For small U, the 
SO bipolaron is approached, and the bandwidth decreases significantly. For reference the tight-binding dispersion normalized 
to the same bandwidth is also plotted. The triplet dispersion is shown in panel (a). 
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FIG. 20: Comparison of the variation of the bandwidth with 
U for the square and triangular lattices. The widening of 
the bandwidth at intermediate U is more significant for the 
square lattice, but is also visible for the triangular lattice. 



In the appendix, we summarize the analytical approach 
to solving the UV-mode\ on a triangular lattice when 
pairs have finite momentum, and the full dispersion can 
be seen in Fig. [2H The energy level crossing is clearly 
visible, and by comparison it is clear that a single s band 



and a single p band are visible in Fig. [2T] Note that there 
are are a total of 6 possible bands associated with the 
bipolaron, although it is unlikely that the higher energy 
bands would be excited in any normal transport process. 

The band crossing can not be seen for bipolarons on the 
square lattice, rather the singlet and triplet pairs become 
degenerate at the zone edge as the phonon frequency in- 
creases. This feature is also visible in the dispersion of 
bipolarons on the chaini^, and is probably a generic fea- 
ture of bipolarons on lattices with cubic symmetry. 



VII. SUMMARY 

In this paper, we computed the properties of bipo- 
larons on two-dimensional lattices using a continuous 
time quantum Monte Carlo algorithm. Properties of 
the bipolaron including the total energy, inverse mass, 
bipolaron radius and number of phonons associated with 
the bipolaron demonstrated the qualitative difference be- 
tween models of electron phonon interaction with long- 
range interaction (screened Frohlich) and those with 
purely local (Holstein) interaction, with only small long 
range tails needed to completely change the properties of 
the bipolaron. A major result of our survey of the pa- 
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FIG. 21: Graphs showing the variation of the singlet and triplet dispersions in the near-neighbor model with u) with A = 8, 
w = not (w = 80t), U = UOt [U = 80t) to match with aU -V model with U = 2At [U = 16t) and V = m {V = 32t). Also 
a) = 30 and w = 6 (and equivalent for the square lattice, a) = 20 and uj = 4). 



rameter space is the existence of extra-light hybrid bipo- 
larons consisting of an on-site and an off-site component 
when on-site and inter-site bipolaron configurations are 
degenerate, similar to the previously reported superlight 
bipolarons, but slightly lighter on the triangular lattice, 
and significantly lighter than the intersite bipolarons on 
the square lattice. We also compute triplet properties 
of the bipolarons. A major surprise is that triplet pairs 
with large momentum (close to the K point) are more 
stable than singlet states on the triangular lattice. 

The potential for the existence of local triplet pairs 

that are preferentially bound has the potential to open 
new avenues in the theory of superconductivity. At low 



momenta, triplet pairs are necessarily higher in energy 
than their singlet counterparts since an exact theorem 
requires the ground state wavefunction of a pair to have 
no nodes. No such theorem exists for states with high 
pair momentum, however the existence of stable triplet 
states is unexpected. Thus, in addition to the poten- 
tial for BCS-BEC crossover as the attractive potential is 
tuned, there is an additional possibility of triplet BCS to 
triplet EEC crossover, or even a transition between sin- 
glet and triplet pairing. We believe that such a possibility 
should be investigated further. 

Our other finding of light hybrid bipolarons on the 
square lattice at intermediate coupling and Coulomb re- 
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pulsion is also significant. While the effect only exists 
over a very narrow range of coupling constants, it shows 
that exotic lattice types are not necessary for small light 
pairs that could have potential to form a high temper- 
ature Bose condensate. As in the case of the triangular 
lattice, light pairs are only a prerequisite for a high tem- 
perature BEC, as other factors such as an absence of 
clustering need to be confirmed. Work to discuss this 
clustering will form the basis of a future publication, but 
clearly there are still some surprises left in the bipolaron 
problem. 
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FIG. 22: The spectrum of bound pair in the UV model on the 
triangular lattice. [/ = 24f and V = 48t, which corresponds 
to A = 8, [/ = 20W and to — > ex. Notice the crossing of the 
lowest singlet and triplet states near the edge of the Brillouin 
zone at the AI point. 



APPENDIX A: UV MODEL ON THE 
TRIANGULAR LATTICE 

The quantum-mechanical problem of two interacting 
particles on a lattice can be solved exactly as long as 
the radius of the interaction is finite. Of special inter- 
est here is the spectrum of bound pairs that form when 
at least part of the potential is attractive. The general 
solution of the bound state problem results in a spectral 
equation2ii2^i2& in the form of a determinant with a size 
equal to the number of lattice sites within the interaction 
range. For the non-retarded UV model with on-site Hub- 
bard repulsion U and nearest-neighbor attraction —V, 
the determinant is 3 x 3 for the one-dimensional chain, 
5x5 for the square lattice, 7 x 7 for the triangular lat- 
tice, and so on. The subsequent analysis can be simplified 
by considering the symmetric and anti-symmetric states 
(that is the singlets and triplets) separately. This effec- 
tively halves the number of interaction sites, and halves 
the size of the determinants as a result. For the UV- 
model on the triangular lattice, the size of the singlet de- 
terminant reduces to 4 (symmetric combinations of the 
nearest neighbors plus the central site), and that of the 
triplet determinant to 3 (anti-symmetric combinations of 
the nearest neighbors). 

Omitting the straightforward but lengthy derivation, 
the spectral equation for the singlet bound pair of the 
lattice momentum K is 
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Here the vectors rij^k denote the nearest neighbor sites 
of the triangular lattice: tiq — (0,0); rii = (1,0); n2 — 

The equation (|Aip can be solved numerically with re- 
spect to the pair energy E for the given pair momentum 
K. For a positive on-site repulsion U and a large enough 
intersite attraction V , the equation has three roots that 
we denote s, di and ^2 states. The evolution of the en- 
ergies as a function of K is shown in figure At the 
high-symmetry points of the Brillouin zone, the system 
of equations (jAip can be further diagonalized into smaller 
blocks by a proper linear transformation. Thus, at the F 
point K = (0, 0), the spectrum splits into the standalone 
ground state s and a degenerate doublet (di, c?2)- In the 
corner of the Brillouin zone, the states s and di become 
degenerate, as can be seen from figure 

The triplet bound states can be treated similarly. The 
exact spectrum equation is a 3 x 3 determinant: 
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with the same meaning of rij.fc and A. The equation does 
not depend on the on-site repulsion at aU, as expected. 
At the r point, the spectrum sphts into a low-energy 
doublet (pi,P2), and a standalone high-energy /-state, 
(see figure [22t . At the iC-point, p2 and / are degenerate. 

Two comments are in order, (i) For strongly- coupled 
pairs, the classification of the different orbital states is 



basically the one of a single particle on a six-site tight- 
binding ring, (ii) Notice the general property of level 
crossing between the lowest singlet and triplet states. 
While the singlet is always the ground state at small 
momenta, the triplet becomes the lowest state at large 
momenta. The same property has been observed for the 
bipolaron, as discussed in the main text of the paper. 
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